Deep convolutional neural network powered terahertz ultra-massive multi-input-multi-output channel estimation method

ABSTRACT

A THz UM-MIMO channel estimation method based on the DCNN comprises the steps: the hybrid spherical and planar-wave modeling (HSPM), by taking a sub-array in the antenna array as a unit, employing the PWM within the sub-array, and employing the SWM among the sub-arrays; estimating the channel parameters between the reference sub-arrays at Tx and Rx through a DCNN, including the angles of departure and arrival, the propagation distance and the path gain; deducing the channel parameters between the reference sub-array and other sub-arrays by utilizing the obtained channel parameters and the geometrical relationships among sub-arrays, and recovering the channel matrix; wherein accurate three-dimensional channel modeling is achieved by the HSPM, which possesses high modeling accuracy and low complexity.

CROSS REFERENCES

This application claims priority to Chinese Patent Application Ser. No. CN202110715991.1 filed on 28 Jun. 2021.

TECHNICAL FIELD

The invention relates to the technical field of THz communications, in particular to the THz UM-MIMO channel modeling and estimation methods based on the DCNN.

BACKGROUND OF THE INVENTION

The accurate three-dimensional channel model is the basis for studying the ultra-massive multi-input-output (UM-MIMO) system, efficiently and accurately estimating the Terahertz (THz) channel information is the premise of realizing the potential of THz systems. However, existing MIMO system research generally employs a planar-wave channel model (PWM). This model is a simplification of the ground-truth spherical-wave channel model (SWM) by approximating signal propagation as a plane. The modeling error can be neglected when the array scale and carrier frequency are small. However, with the increment of the dimension of the antenna array and the carrier frequency, the accuracy of the planar-wave model is greatly reduced. Therefore, it is necessary to consider the most accurate SWM. However, the complexity of the SWM is positively related to the number of antennas. As the number of antennas in the THz UM-MIMO system increases, the complexity of the SWM increases correspondingly.

Existing channel estimation (CE) methods of the MIMO systems usually adopt a PWM, which is a simplification of the SWM. In the THz UM-MIMO systems, due to the increment of the dimension of the antenna array and the carrier frequency, the accuracy of the PWM is greatly reduced, and the spherical-wave should be considered. However, with the change of the channel model, existing CE algorithms lose their effectiveness.

SUMMARY OF THE INVENTION

Aiming at the problem that the channel model of the existing THz UM-MIMO system is inaccurate and the CE method is to be studied, the invention provides a THz UM-MIMO CE method based on a deep convolutional neural network (DCNN). The method exploits a hybrid spherical and planar-wave modeling (HSPM) for accurate three-dimensional channel modeling, which realizes high modeling accuracy and low complexity, and proposes a matched CE method according to the established HSPM, to effectively obtain the channel information and achieve optimal resource allocation, and fully exert the potential of the THz communication systems.

The present invention is implemented by the following technical solutions:

The invention relates to a THz UM-MIMO CE method based on a HSPM, comprising the following steps:

Step 1: the HSPM channel modeling, within a sub-array, the PWM is adopted, and the SWM is adopted among the sub-arrays.

Step 2: estimating the angles of departure and arrival, the propagation distance and the path gain between the reference sub-arrays in the step 1 by the DCNN.

Step 3: The channel parameters between the reference sub-array and the remaining sub-arrays are deduced by using the channel parameters obtained in step 2 and the geometric relationships among the sub-arrays, and the channel matrix is recovered.

Technical Effect

The present invention solves the problem in the prior art that the PWM of the UM-MIMO system is not accurate enough, the complexity of the SWM is too high, and the CE algorithm matched with the UM-MIMO system is not effective enough.

Compared with the prior art, the HSPM is adopted, the PWM is adopted in the sub-array, a SWM is adopted among different sub-arrays, and the high precision and low complexity channel modeling is realized. At the same time, through the two-stage CE mechanism, the channel parameter estimation of the reference sub-array is completed by designing the DCNN, and then the residual parameter estimation is completed by using the geometric relationship of the channel parameters between the reference sub-arrays, and the channel matrix is recovered. Since only a DCNN needs to be used for parameter estimation of the reference sub-array, the required number of channel parameters to be estimated is the same as that by using the PWM. Since the parameter estimation at the second stage can be performed in parallel, the proposed method achieves low complexity and high effectiveness of CE compared to conventional CE algorithms

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of a system according to the present invention;

FIG. 2 is a structural diagram of the DCNN according to the present invention;

FIG. 3A is a schematic diagram of the effect of an embodiment: the error of different channel models at different communication distances is exploited;

FIG. 3B is a schematic diagram of the effect of an embodiment: the error of different channel models at different sub-array spacing is determined

FIG. 3C is a schematic diagram of the effect of an embodiment: the error of different channel models at different carrier frequencies;

FIG. 4 is a schematic performance comparison diagram of the method of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

As shown in FIG. 1 , a THz UM-MIMO CE system based on a DCNN according to an embodiment includes: a radio frequency (RF)-chain, an analog beamformer, an analog combiner, and a digital combiner. The RF-chain at the transmitter (Tx) transmits the baseband pilot signal and perform digital beamforming, to obtain the RF signal. The analog beamformer performs analog beamforming according to a predefined beam codebook, and transmits the analog signal to the channel. The analog combiner at the receiver (Rx) combines the received signals and outputs it to the digital combiner. The digital combiner performs signal processing on the combined analog signal to obtain a baseband signal. Finally, the baseband processes the digital signal to complete the CE.

The present embodiment relates to a THz UM-MIMO CE method based on a DCNN is composed of the following steps:

Step 1, as shown in FIG. 1 , take a sub-array as a unit, the PWM is adopted in the sub-array, and the SWM is used between the sub-arrays to complete the HSPM, specifically it is composed of:

-   -   1) dividing the antennas at Tx and Rx into K _(t) and K _(r)         sub-arrays, respectively, and different sub-arrays have the same         number of multi-path N _(p). The channel gains between different         sub-arrays have the same amplitude, while the phase of the         channel gain is changed due to different geometric distances and         transceiver angles. The HSPM can be expressed as

${H_{HSPM} = {\sum_{p = 1}^{N_{p}}\begin{bmatrix} {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{11}}{a_{rp}^{11}\left( a_{tp}^{11} \right)}^{H}} & \ldots & {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{1K_{t}}}{a_{rp}^{1K_{t}}\left( a_{tp}^{1K_{t}} \right)}^{H}} \\ \ldots & & \ldots \\ {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{K_{r}1}}{a_{rp}^{K_{r}1}\left( a_{tp}^{K_{r}1} \right)}^{H}} & \ldots & {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{K_{r}1}}{a_{rp}^{K_{r}K_{t}}\left( a_{tp}^{K_{r}K_{t}} \right)}^{H}} \end{bmatrix}}},$

where N _(p) is the number of multi-paths in the channel, α_(p) ¹¹ is the channel gain of the p^(th) path, λ is the signal wavelength, D_(p) ^(k) ^(r) ^(k) ^(t) is the transmission distance of the p^(th) path, α_(rp) ^(k) ^(r) ^(k) ^(t) and α_(tp) ^(k) ^(r) ^(k) ^(t) are the array steering vectors at Tx and Rx, respectively, whose value are determined by the propagation angles.

The block structure of the HSPM H _(HSPM) refers to: each block

${❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{k_{r}k_{t}}}{a_{rp}^{k_{r}k_{t}}\left( a_{tp}^{k_{r}k_{t}} \right)}^{H}$

is a PWM matrix between the sub-arrays at Tx and Rx. Among between different sub-arrays, the SWM is considered, resulting in different transceiver angles, and thus α_(rp) ^(k) ^(r) ^(k) ^(t) and α_(tp) ^(k) ^(r) ^(k) ^(t) are different.

-   -   2) According to the channel model, the received signal can be         expressed as Y=W ^(H) H _(HSPM) F+N, W is the received combining         matrix comprising the analog combining and digital combining         matrices, and F is a transmitted beamforming matrix comprising         an analog beamforming matrix and a digital beamforming matrix,         and N is a noise signal.

Step 2, using the first sub-array at Tx and Rx as the reference sub-arrays, and perform parameter estimation by using the DCNN, including the departure angles, the angles of arrival the propagation distance, and the path gain.

2.1 The DCNN estimates the departure angles (θ_(t, p), ϕ_(t, p)), the angle of arrival (θ_(r, p), ϕ_(r, p)), the propagation D _(p) ^(k) ^(r) ^(k) ^(t) , and the path gain |α_(p) ¹¹| according to Re{Y}, Re{Y} and |Y|.

As shown in FIG. 2 , there are totally 15 layers in the DCNN, which includes an input layer, seven convolutional layers (CV), four max-pooling layers (MP), a flattening layer, and a fully connected (FC) output layer, wherein the input layer receives the element real value Re{Y} of the channel observation matrix Y, the element imaginary value lm{Y} and the element absolute value |Y|; the CV respectively comprise 16, 32, 64, 128, 62, 32, 16 convolution kernels, wherein a convolution filter with a size of 3×3 is deployed to extract features; zero padding (ZP) and batch normalization (BN) are deployed in the convolution process; the pooling layer uses the maximum pooling method to extract the maximum value in the 2×2 size pool to reduce the number of network dimensions, and simplify the training process. The flattening layer rearranges the neurons into one dimension and is connected to the FC output layer for outputting the estimated emission angle, the angle of arrival, the propagation distance, and the path gain.

Step 3, deriving channel parameters between the reference sub-array and the remaining sub-arrays by using the geometric relationship between the sub-arrays obtained in Step 2, and restoring the channel matrix, specifically comprising:

1) Calculating parameters of the line-of-sight path: projecting the line-of-sight path as shown in FIG. 1 to the X-Y plane and the Y-Z plane, respectively, and deriving the

${{{angles}\theta_{t}^{k_{t}k_{r}}} = {\arccos\left( \frac{D_{xy}^{11}\cos\theta_{t}^{11}}{D_{xy}^{k_{t}k_{r}}} \right)}},{\theta_{r}^{k_{t}k_{r}} = {\theta_{r}^{11} + {\arcsin\left( \frac{\Delta d_{x}\cos\theta_{t}^{11}}{D_{xy}^{k_{t}k_{r}}} \right)}}},$ ${\phi_{t}^{k_{t}k_{r}} = {\arccos\left( \frac{D_{yz}^{11}\cos\phi_{t}^{11}}{D_{yz}^{k_{t}k_{r}}} \right)}},{\phi_{r}^{k_{t}k_{r}} = {{- \phi_{r}^{11}} + {\arcsin\left( \frac{\Delta d_{z}\cos\theta_{t}^{11}}{D_{yz}^{k_{t}k_{r}}} \right)}}}$

and the propagation distance

${D^{k_{t}k_{r}} = \frac{D_{yz}^{k_{t}k_{r}}}{\cos\theta_{t}^{11}}},$

wherein: θ_(t) ^(k) ^(t) ^(k) ^(r) and ϕ_(t) ^(k) ^(t) ^(k) ^(r) denote the received azimuth and elevation angles of the k_(t) sub-array at Tx end and the k_(r) sub-array at Rx, respectively.

${D_{xy}^{11} = {{D^{11}\cos\phi_{t}^{11}D_{yz}^{11}} = {D^{11}\cos\theta_{t}^{11}}}},{D_{xy}^{k_{t}k_{r}} = \sqrt{\left( {\Delta d_{x}^{k_{t}k_{r}}} \right)^{2} + \left( D_{xy}^{11} \right)^{2} - {2\Delta d_{x}^{k_{t}k_{r}}D_{xy}^{11}\sin\theta_{t1}^{11}}}},$ ${D_{yz}^{k_{t}k_{r}} = \sqrt{\left( {\Delta d_{z}^{k_{t}k_{r}}} \right)^{2} + \left( D_{yz}^{11} \right)^{2} - {2\Delta d_{z}^{k_{t}k_{r}}D_{yz}^{11}\sin\phi_{r}^{11}}}},$

where D ¹¹ represent the distance between the reference sub-arrays at Tx and Rx. Δd_(x) ^(k) ^(t) ^(k) ^(r) refer to the relative displacement of the k_(r) sub-array at Rx to the k_(t) sub-array at Tx along the X-axis. Δd_(z) ^(k) ^(t) ^(k) ^(r) is the relative displacement along the Z-axis of the k_(r) sub-array at Rx to the k_(t) sub-array at Tx.

${D_{yz}^{11} = {D^{11}\cos\theta_{t}^{11}}},{D_{xy}^{k_{t}k_{r}} = {\sqrt{\left( {\Delta d_{x}^{k_{t}k_{r}}} \right)^{2} + \left( D_{xy}^{11} \right)^{2} - {2\Delta d_{x}^{k_{t}k_{r}}D_{xy}^{11}\sin\theta_{t1}^{11}}}.}}$

2) calculating the plane equation of the reflecting surface: solving the equation A_(p)x+B_(p)y+C_(p)z+D_(p)=0 of the reflecting surface by means of the coordinates of the transceiver end and the equation of the equation,

${A_{p} = \frac{\sin{\phi_{tp}^{11}\left( {{C_{p}^{11}\cos\theta_{sp}^{11}\cos\phi_{sp}^{11}} - {B_{p}^{11}\sin\phi_{tp}^{11}}} \right)}}{\cos\phi_{sp}^{11}\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\cos\theta_{sp}^{11}} - {B_{p}^{11}\sin\theta_{p}^{11}}} \right)}}},{B_{p} = \frac{\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\sin\phi_{sp}^{11}} - {C_{p}^{11}\sin\theta_{sp}^{11}\cos\phi_{sp}^{11}}} \right)}}{\cos\phi_{sp}^{11}\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\cos\theta_{sp}^{11}} - {B_{p}^{11}\sin\theta_{sp}^{11}}} \right)}}},$

C_(p)=1, wherein the lower corner mark _(p) indexes the propagation path, and the remaining parameter definitions are the same as the before.

3) obtaining parameters of the non-line-of-sight path after obtaining the plane equation of the reflecting surface, specifically comprising:

${\theta_{tp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{x_{p}^{k_{r}k_{r}}}{\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2}}} \right\rbrack}},{\phi_{tp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{z_{p}^{k_{r}k_{r}}}{\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2} + \left( z_{p}^{k_{r}k_{r}} \right)^{2}}} \right\rbrack}},$ ${\theta_{rp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{x_{p}^{k_{r}k_{r}} - R_{k_{r}x}}{\sqrt{\left( {x_{p}^{k_{r}k_{r}} - R_{1x}} \right)^{2} + \left( {y_{p}^{k_{r}k_{r}} - R_{1y}} \right)^{2}}} \right\rbrack}},$ $\phi_{rp}^{k_{r}k_{r}} = {{\arcsin\left\lbrack \frac{z_{p}^{k_{r}k_{r}} - R_{k_{r}z}}{\sqrt{\left( {x_{p}^{k_{r}k_{r}} - R_{k_{r}x}} \right)^{2} + \left( {y_{p}^{k_{r}k_{r}} - R_{k_{r}y}} \right)^{2} + \left( {z_{p}^{k_{r}k_{r}} - R_{k_{r}z}} \right)^{2}}} \right\rbrack}.}$ where ${D_{p}^{k_{t}k_{r}} = {\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2} + \left( z_{p}^{k_{r}k_{r}} \right)^{2}} + \sqrt{\left( {R_{k_{r}x} - x_{p}^{k_{r}k_{r}}} \right)^{2} + \left( {R_{k_{r}y} - y_{p}^{k_{r}k_{r}}} \right)^{2} + \left( {R_{k_{r}z} - z_{p}^{k_{r}k_{r}}} \right)^{2}}}},x_{p}^{k_{r}k_{r}},y_{p}^{k_{r}k_{r}},z_{p}^{k_{r}k_{r}}$

denotes the coordinates of the remaining subarrays od the reflection plane, (R_(k) _(r) _(x), R_(k) _(r) _(y), R_(k) _(r) _(z))=(D¹¹ sin θ_(t1) ¹¹ cos ϕ_(t1) ¹¹+d_(rx) ^(k) ^(r) , D¹¹ cos θ_(t1) ¹¹ cos ϕ_(t1) ¹¹, D¹¹ sin ϕ_(t1) ¹¹−d_(rz) ^(k) ^(r) ).

In a specific practical experiment, we set 1024 antennas and 4 subarrays at Rx, the communication distance is 20 m and the carrier frequency is 0.3 THz, the precision of the HSPM combined channel model is improved by 14 dB. As shown in FIG. 3A-C and FIG. 4 and Table 1, compared to the existing CE method, the CE method based on the DCNN can complete CE in 0.172 ms, and achieves an estimation precision of 5.2 dB.

TABLE 1 Method Computational Complexity Running Time (ms) OMP [16]

 ((N_(p)N_(t))²) 221 AMP [17]

 ((N_(p)N_(t))²) 372 CNN 07

 (b(N_(t)N_(r))²) 3.64 RNN [40]

 (cC³KtK_(r)) 0.085 DCNN

 (b(C²KtK_(r))² 0.172

Compared with the prior art, the HSPM provided by the invention achieves high precision and low complexity of channel modeling, and compared with the PWM, the precision is improved by 14 dB, and compared with the SWM, the complexity is reduced by 99%. The proposed CE method based on the DCNN achieves the improvement of the CE precision of 5.2 dB, and at the same time, since only the channel parameters between the reference sub-arrays need to be estimated, the CE overhead of the method is reduced by 93%.

The above specific implementation may be locally adjusted by a person skilled in the art without departing from the principle and spirit of the present disclosure, and the protection scope of the present disclosure is subject to the claims and is not limited by the specific embodiments described above, and various implementations within the scope of the present disclosure are not limited by the present disclosure. 

What is claimed is:
 1. A THz UM-MIMO CE method based on a DCNN, comprising the following steps: step i), the HSPM modeling, which takes the sub-array as a unit, using the PWM in the sub-array, and models the channel among sub-arrays by the SWM; step ii), using the first sub-array at the transceiver end as a reference sub-array, using a DCNN to estimate the departure angle, the angle of arrival, the propagation distance and the path gain between the reference sub-arrays according to real values, element imaginary values and element absolute values of the channel observation matrix; step iii), deriving channel parameters between the reference sub-array and the remaining sub-arrays by using the channel parameters obtained in step ii) and the geometric relationship between the sub-arrays, and reconstruct the channel matrix.
 2. The THz UM-MIMO CE method according to claim 1, wherein step i) comprises the following steps: a) dividing the antennas at Tx and Rx into K _(t) and K _(r) sub-arrays, respectively, and different sub-arrays have the same multi-path number A the amplitude of the channel gain between different sub-arrays is the same, while the phase of the channel gain is changed due to different geometric distances and transceiver angles, to obtain the block structured channel model: ${H_{HSPM} = {\overset{N_{p}}{\sum\limits_{p = 1}}\begin{bmatrix} {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{11}}{a_{rp}^{11}\left( a_{tp}^{11} \right)}^{H}} & \ldots & {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{1K_{t}}}{a_{rp}^{1K_{t}}\left( a_{tp}^{1K_{t}} \right)}^{H}} \\ \ldots & & \ldots \\ {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{K_{r}1}}{a_{rp}^{K_{r}1}\left( a_{tp}^{K_{r}1} \right)}^{H}} & \ldots & {{❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{K_{r}1}}{a_{rp}^{K_{r}K_{t}}\left( a_{tp}^{K_{r}K_{t}} \right)}^{H}} \end{bmatrix}}},$ wherein: N_(p) is the number of multi-path in the channel, α_(p) ¹¹ is the channel gain of the p-th path, _(λ) is the signal wavelength, D_(p) ^(k) ^(r) ^(k) ^(t) is the transmission distance of the p-th path, a_(rp) ^(k) ^(r) ^(k) ^(t) and a_(tp) ^(k) ^(r) ^(k) ^(t) the antenna array guide vector for the transceiver end, respectively; and b) receiving a signal according to a channel model, wherein the received signal Y=W ^(H)H_(HSPM) F+N, W is a received combining matrix comprising an analog combining matrix and a digital combining matrix, and F is a transmitting beamforming matrix comprising an analog beamforming matrix and a digital beamforming matrix, and N is a noise signal.
 3. The THz UM-MIMO CE method according to claim 2, wherein the block structured channel model refers to: each block ${❘\alpha_{p}^{11}❘}e^{{- j}\frac{2\pi}{\lambda}D_{p}^{k_{r}k_{t}}}{a_{rp}^{k_{r}k_{t}}\left( a_{tp}^{k_{r}k_{t}} \right)}^{H}$ is a channel matrix between sub-arrays at the Tx and Rx, the channel matrix is the PWM matrix within the block, and among different sub-arrays, the channel is spherical-wave modeled, resulting in different transceiver angles, and thus the a_(rp) ^(k) ^(r) ^(k) ^(t) and a_(tp) ^(k) ^(r) ^(k) ^(t) are different.
 4. The THz UM-MIMO CE method according to claim 1, wherein the DCNN comprises an input layer, seven convolutional layers (CV), four max-pooling layers (MP), a flattening layer, and a fully connected (FC) output layer, wherein the input layer receives the element real value Re{Y} of the channel observation matrix Y, the element imaginary value lm{Y} and the element absolute value |Y|; the CV respectively comprise 16, 32, 64, 128, 62, 32, 16 convolution kernels, wherein a convolution filter with a size of 3×3 is deployed to extract features; zero padding (ZP) and batch normalization (BN) are deployed in the convolution process; the pooling layer uses the maximum pooling method to extract the maximum value in the 2×2 size pool to reduce the number of network dimensions, and simplify the training process. The flattening layer rearranges the neurons into one dimension and is connected to the FC output layer for outputting the estimated emission angle, the angle of arrival, the propagation distance, and the path gain.
 5. The THz UM-MIMO CE method according to claim 1, wherein the step iii) comprises: a) calculating parameters of the line-of-sight path: respectively projecting the line-of-sight path to the X-Y plane and the Y-Z plane, and deriving a transceiver angle between the remaining sub-arrays under the line of sight ${\theta_{t}^{k_{t}k_{r}} = {\arccos\left( \frac{D_{xy}^{11}\cos\theta_{t}^{11}}{D_{xy}^{k_{r}k_{t}}} \right)}},{\theta_{r}^{k_{t}k_{r}} = {\theta_{r}^{11} + {\arcsin\left( \frac{\Delta d_{x}\cos\theta_{t}^{11}}{D_{xy}^{k_{t}k_{r}}} \right)}}},$ ${\phi_{t}^{k_{t}k_{r}} = {\arccos\left( \frac{D_{yz}^{11}\cos\phi_{t}^{11}}{D_{yz}^{k_{t}k_{r}}} \right)}},{\phi_{r}^{k_{t}k_{r}} = {{- \phi_{r}^{11}} + {\arcsin\left( \frac{\Delta d_{z}\cos\theta_{t}^{11}}{D_{yz}^{k_{t}k_{r}}} \right)}}}$ and a propagation distance ${D^{k_{t}k_{r}} = \frac{D_{yz}^{k_{t}k_{r}}}{\cos\theta_{t}^{11}}},$ wherein: θ_(t) ^(k) ^(t) ^(k) ^(r) and ϕ_(t) ^(k) ^(t) ^(k) ^(r) are the receiving azimuth angle and the elevation angle of the k_(r) sub-array at the k_(t) transmitting end, while θ_(r) ^(k) ^(t) ^(k) ^(r) and ϕ_(r) ^(k) ^(t) ^(k) ^(r) the k_(r) sub-array at the k_(t) receiving end, respectively, ${D_{xy}^{11} = {{D^{11}\cos\phi_{t}^{11}D_{yz}^{11}} = {D^{11}\cos\theta_{t}^{11}}}},{D_{xy}^{k_{t}k_{r}} = \sqrt{\left( {\Delta d_{x}^{k_{t}k_{r}}} \right)^{2} + \left( D_{xy}^{11} \right)^{2} - {2\Delta d_{x}^{k_{t}k_{r}}D_{xy}^{11}\sin\theta_{t1}^{11}}}},{D_{yz}^{k_{t}k_{r}} = \sqrt{\left( {\Delta d_{z}^{k_{t}k_{r}}} \right)^{2} + \left( D_{yz}^{11} \right)^{2} - {2\Delta d_{z}^{k_{t}k_{r}}D_{yz}^{11}\sin\phi_{r}^{11}}}},$ wherein: D¹¹ is the distance between the transmitting end and the reference sub-array of the receiving end, and Δd_(x) ^(k) ^(t) ^(k) ^(r) the relative displacement along the X axis of the k_(r) sub-array at the k_(t) transmitting end, here D ¹¹ represent the distance between the reference sub-arrays at Tx and Rx. Δd_(x) ^(k) ^(t) ^(k) ^(r) refer to the relative displacement of the k_(r) sub-array at Rx to the k_(t) sub-array at Tx along the X-axis. Δd_(z) ^(k) ^(t) ^(k) ^(r) is the relative displacement along the Z-axis of the k_(r) sub-array at Rx to the k_(t) sub-array at Tx; D _(yz) ¹¹ =D ¹¹ cos θ_(t) ¹¹, ${D_{xy}^{k_{t}k_{r}} = \sqrt{\left( {\Delta d_{x}^{k_{t}k_{r}}} \right)^{2} + \left( D_{xy}^{11} \right)^{2} - {2\Delta d_{x}^{k_{t}k_{r}}D_{xy}^{11}\sin\theta_{t1}^{11}}}};$ b) calculating the plane equation of the reflecting surface: solving the equation A_(p)x+B_(p)y+C_(p)z+D_(p)=0, by means of the reflecting surfaces by means of the coordinates of the transceiver end and the law of reflection, wherein: each parameter is: ${A_{p} = \frac{\sin{\phi_{tp}^{11}\left( {{C_{p}^{11}\cos\theta_{sp}^{11}\cos\phi_{sp}^{11}} - {B_{p}^{11}\sin\phi_{tp}^{11}}} \right)}}{\cos\phi_{sp}^{11}\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\cos\theta_{sp}^{11}} - {B_{p}^{11}\sin\theta_{p}^{11}}} \right)}}},{B_{p} = \frac{\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\sin\phi_{sp}^{11}} - {C_{p}^{11}\sin\theta_{sp}^{11}\cos\phi_{sp}^{11}}} \right)}}{\cos\phi_{sp}^{11}\sin{\phi_{tp}^{11}\left( {{A_{p}^{11}\cos\theta_{sp}^{11}} - {B_{p}^{11}\sin\theta_{sp}^{11}}} \right)}}},{C_{p} = 1},$ the lower corner mark _(p) index the propagation path; and c) obtaining parameters of the non-line-of-sight path after obtaining the plane equation of the reflecting surface, ${\theta_{tp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{x_{p}^{k_{r}k_{r}}}{\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2}}} \right\rbrack}},{\phi_{tp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{z_{p}^{k_{r}k_{r}}}{\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2} + \left( z_{p}^{k_{r}k_{r}} \right)^{2}}} \right\rbrack}},{\theta_{rp}^{k_{t}k_{r}} = {\arcsin\left\lbrack \frac{x_{p}^{k_{r}k_{r}} - R_{k_{r}x}}{\sqrt{\left( {x_{p}^{k_{r}k_{r}} - R_{1x}} \right)^{2} + \left( {y_{p}^{k_{r}k_{r}} - R_{1y}} \right)^{2}}} \right\rbrack}},{\phi_{rp}^{k_{r}k_{r}} = {{\arcsin\left\lbrack \frac{z_{p}^{k_{r}k_{r}} - R_{k_{r}z}}{\sqrt{\left( {x_{p}^{k_{r}k_{r}} - R_{k_{r}x}} \right)^{2} + \left( {y_{p}^{k_{r}k_{r}} - R_{k_{r}y}} \right)^{2} + \left( {z_{p}^{k_{r}k_{r}} - R_{k_{r}z}} \right)^{2}}} \right\rbrack}.}}$ wherein ${D_{p}^{k_{t}k_{r}} = {\sqrt{\left( x_{p}^{k_{r}k_{r}} \right)^{2} + \left( y_{p}^{k_{r}k_{r}} \right)^{2} + \left( z_{p}^{k_{r}k_{r}} \right)^{2}} + \sqrt{\left( {R_{k_{r}x} - x_{p}^{k_{r}k_{r}}} \right)^{2} + \left( {R_{k_{r}y} - y_{p}^{k_{r}k_{r}}} \right)^{2} + \left( {R_{k_{r}z} - z_{p}^{k_{r}k_{r}}} \right)^{2}}}},x_{p}^{k_{r}k_{r}},y_{p}^{k_{r}k_{r}},$ z_(p) ^(k) ^(r) ^(k) ^(r) are the coordinates of the reflection point of the remaining sub-array on the reflecting surface (R_(k) _(r) _(x), R_(k) _(r) _(y), R_(k) _(r) _(z))=(D ¹¹ sin θ_(t1) ¹¹ cos ϕ_(t1) ¹¹+d_(rx) ^(k) ^(r) , D ¹¹ cos θ_(t1) ¹¹ cos ϕ_(t1) ¹¹ , D ¹¹ sin ϕ_(t1) ¹¹−d_(rz) ^(k) ^(r) ), and d_(rx) ^(k) ^(r) is the distance from k_(r) th sub-array of the receiving end to the reference sub-array along the X-axis, d_(rz) ^(k) ^(r) is the distance from the k_(r) th sub-array of the receiving end to the reference sub-array along the Z-axis.
 6. A THz UM-MIMO CE system realizing the said THz UM-MIMO CE method according to either claim 1 comprising: a RF-chain, an analog beamformer, an analog combiner, and a digital combiner, wherein the RF-chain of the transmitting end receives the baseband pilot signal to perform digital beamforming to obtain a radio frequency signal, and the analog beamformer performs analog beamforming according to a preset beam codebook signal and transmits the analog beamforming signal to the channel; and the analog combiner at Rx end combines the received signals and then outputs the received signals to the digital beam forming unit; the digital beamformer performs beamforming processing according to the combined analog signal to obtain a digital baseband signal, and the baseband end processes the digital baseband signal to complete the CE. 